setwd("E:/5hmc_file/2_5hmc_yjp_bam/ASM/20201216")
library(ggplot2)
file=read.csv("4.Coordinates.site.csv",header=T)
id=file$unitID
setwd("E:/1.甲基化分析/ASM/ASM_snp-onlyWGS")
sel1=list.files(pattern = ".snp")

i=2
file=read.table(sel1[i],head=T,sep="\t")
file=file[file$exits=="T",]
file$unitID=paste(file$chrom,file$position,file$ref,file$var,sep=":")
file=file[file$unitID %in% id,]
if(dim(file)[1]>0){
  file$tumor_var_freq=as.numeric(gsub("%","",file$tumor_var_freq))/100
  file$normal_var_freq=as.numeric(gsub("%","",file$normal_var_freq))/100
  file=data.frame(file$unitID,file$normal_var_freq,file$tumor_var_freq)
  str=unlist(strsplit(gsub(".snp","",sel1[i]),"_"))
  names(file)=c("unitID",paste0(str[1],".var.freq"),paste0(str[2],".var.freq"))
  rt=file
}


for(i in 3:length(sel1)){
  file=read.table(sel1[i],head=T,sep="\t")
  file=file[file$exits=="T",]
  file$unitID=paste(file$chrom,file$position,file$ref,file$var,sep=":")
  file=file[file$unitID %in% id,]
  if(dim(file)[1]>0){
    file$tumor_var_freq=as.numeric(gsub("%","",file$tumor_var_freq))/100
    file$normal_var_freq=as.numeric(gsub("%","",file$normal_var_freq))/100
    file=data.frame(file$unitID,file$normal_var_freq,file$tumor_var_freq)
    str=unlist(strsplit(gsub(".snp","",sel1[i]),"_"))
    names(file)=c("unitID",paste0(str[1],".var.freq"),paste0(str[2],".var.freq"))
    rt=merge(rt,file,by="unitID",all = T)
  }
}


setwd("E:/5hmc_file/2_5hmc_yjp_bam/ASM/20201216")
library(ggplot2)
file=read.csv("4.Coordinates.site.csv",header=T)
filea=read.csv("../20201112做汇总表/all.FDR.sig.at.least.one.add.direction.same.diff.csv",head=T)
ID2=c("X1T","X2B","M7","M8","M5","M6","M1","M2","M41","M42","M43","M44","M47","M48","M3","M4","M11","M12","M27","M28","M29","M30","M25","M26","M36","M35","M32","M31","M17","M18","M19","M20","M21","M22","M39","M40","M49","M50","M51","M52")
ID=paste(ID2[seq(2,40,2)],ID2[seq(1,40,2)],sep="_")
status=c(rep("DC",9),rep("CC",5),rep("HC",4),rep("DC",2))
sample_info=data.frame(ID,status)
row.names(sample_info)=sample_info$ID
file=rt
allsample=colnames(file)[grep("var.freq",colnames(file))]###AShM
sel1=allsample

for(i in 1:2){
mta1=as.numeric(file[i,c(sel1[seq(2,length(sel1),2)])])
mta2=as.numeric(file[i,c(sel1[seq(1,length(sel1),2)])])
sample_name=gsub(".var.freq","",colnames(file)[grep(".var.freq",colnames(file))])
sample_name=paste(sample_name[seq(1,16,2)],sample_name[seq(2,16,2)],sep="_")
sample_name=sample_name[!is.na(mta1)]

mta1=mta1[!is.na(mta1)]
mta2=mta2[!is.na(mta2)]

group=as.character(sample_info[sample_name,]$status)
df_con=data.frame(VAF=mta2,type="control",sample=sample_name,group)
df_cas=data.frame(VAF=mta1,type="case",sample=sample_name,group)
df=rbind(df_con,df_cas)
df=df[df$group=="DC",]

j=file[i,]$unitID
fileMF=paste(filea[filea$unitID==j,]$avsnp150,"_",filea[filea$unitID==j,]$Gene.refGene,".ASM.pdf",sep="")

p=ggplot(data=df,aes(x=type,y=VAF,group=sample,color=group))+
  geom_line(size=1)+geom_point(aes(shape=type),size=3)+
  ylim(0,1)+theme_light()+labs(x="ASM",y="variant allele frequency",title = fileMF)+
  theme(axis.title.y=element_text(size=16),axis.text.x=element_text(face="bold",size=12),axis.title.x=element_text(size=16))

pdf(fileMF,width = 5,height = 5)
print(p)
dev.off()
}
